home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1997 August / Walnut Creek CDROM.7z / VOL_400 / 446_01 / DOC / FDMSTART / EX7 / CONSLAW.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-04-18  |  6.0 KB  |  199 lines
  1. #include <ConsLaw.h>
  2. #include <MenuSystem.h>
  3. #include <createHCLSchemes.h>
  4. #include <createFluxFunc.h>
  5. #include <createInitFunc.h>
  6. #include <DrawFD.h>
  7.  
  8. ConsLaw:: ConsLaw ()
  9. {  p_value = 0;  /* gives p=1 field due to exp(0) */ }
  10.  
  11. BooLean ConsLaw:: ok () const
  12. {
  13.   BooLean b = dpTRUE;
  14.   HANDLE_OK("ConsLaw::ok",u,b)
  15.   HANDLE_OK("ConsLaw::ok",u_prev,b)
  16.   HANDLE_OK("ConsLaw::ok",p,b)
  17.   if (!scheme.ok()) { b = dpFALSE;
  18.     errorFP("ConsLaw::ok","scheme is not rebinded to a field (empty handle)");}
  19.   if (!flux.ok()) { b = dpFALSE;
  20.     errorFP("ConsLaw::ok","flux is not rebinded (empty handle)"); }
  21.   if (time_points_for_plot.getNoMembers() < 1)
  22.     warningFP("ConsLaw::ok","no plots will be made! (no time points given)");
  23.   return b;
  24. }
  25.  
  26. void ConsLaw:: adm (MenuSystem& menu)
  27. {
  28.   define (global_menu);   // define menu items
  29.   global_menu.prompt();   // prompt the user (read input data)
  30.   scan (global_menu);     // load user's input data into the class variables
  31. }
  32.  
  33. void ConsLaw:: solveProblem ()
  34. {
  35.   if (!ok())
  36.     fatalerrorFP("ConsLaw::driver",
  37.          "object is not properly initialized, have you called adm?");
  38.   setIC ();            // set initial conditions
  39.   timeLoop ();         // perform time stepping
  40. }
  41.  
  42. void ConsLaw:: define (MenuSystem& menu, int level)
  43. {
  44.   menu.addItem (level,       // menu level (main, submenu, subsubmenu)
  45.                 "grid",      // item name, used as command
  46.                 "grid",      // corresponding Unix cml arg., +grid
  47.                 "GridLattice::scan syntax",  // help/description
  48.                 "d=1, domain: [0,1] index: [0:20]",  // default answer
  49.                 "S",         // indicates that valid answer is a string
  50.                 'g', 'g');   // short forms (-g instead of +grid etc)
  51.   menu.addItem (level,
  52.                 "time parameters",
  53.                 "timeprm",
  54.                 "TimePrm::scan syntax",
  55.                 "dt=0.05, t in [0,0.9]",
  56.                 "S",
  57.                 't','t');
  58.   menu.addItem (level,
  59.                 "time points for plot",
  60.                 "plot_times",
  61.                 "Set::scan syntax, eoi: <CR>",
  62.                 "0.2 0.5 0.8",
  63.                 "S", // improvement: could be intelligent from t in [0:?]
  64.                 'p','p');
  65.   menu.addItem (level,
  66.                 "scheme",
  67.                 "scheme",
  68.                 "classname in HCLSchemes hierarchy",
  69.                 "UpwindScheme1D",
  70.                 validationString(hierHCLSchemes()),
  71.                 's','s');
  72.  
  73.   menu.addItem (level,
  74.                 "flux function",
  75.                 "flux",
  76.                 "functor, name in FluxFunc hierarchy",
  77.                 "FluxLinear",
  78.                 validationString(hierFluxFunc()),
  79.                 'f','f');
  80.   menu.addItem (level,
  81.                 "initial function",
  82.                 "u0",
  83.                 "functor, name in InitFunc hierarchy",
  84.                 "InitConst1",
  85.                 validationString(hierInitFunc()),
  86.                 '0','0');
  87.   menu.addItem (level,
  88.                 "uL",
  89.                 "uL",
  90.                 "value of u(x,t) at the left boundary",
  91.                 "1",
  92.                 "R1",
  93.                 'l','l');
  94.   menu.addItem (level,
  95.         "p-mean",
  96.         "E[p]",
  97.         "mean value of the p-field (f(u;p))",
  98.         "0",
  99.         "R1",
  100.         'P','P');
  101. }
  102.  
  103. void ConsLaw:: scan (MenuSystem& menu)
  104. {
  105.   grid.rebind (new GridLattice(1));  grid->scan (menu.get("grid"));
  106.   u.rebind (new FieldFD (grid(), "u"));   
  107.   u_prev.rebind (new FieldFD (grid(), "u_prev"));
  108.   tip.scan (menu.get("time parameters"));
  109.   // add the end mark ';' to the end (required by TimePrm::scan):
  110.   time_points_for_plot.scan (menu.get("time points for plot") + " ;");
  111.   flux_tp = menu.get("flux function");  
  112.   flux.rebind (createFluxFunc(flux_tp,*this));
  113.  
  114.   u0_tp = menu.get("initial function");
  115.   u0.rebind (createInitFunc(u0_tp,*this));
  116.  
  117.   uL = menu.get("uL").getReal();
  118.   u->valueIndex(0) = u_prev->valueIndex(0) = uL;
  119.  
  120.   scheme_tp = menu.get("scheme");  // scheme type
  121.   scheme.rebind (createHCLSchemes (scheme_tp, *this));
  122.  
  123.   // p-field, set p = exp(p_value):
  124.   p.rebind (new FieldFD (grid(), "p")); 
  125.   p_value = menu.get("p-mean").getReal();
  126.   p->fill (exp(p_value));
  127.  
  128.   file.open (CaseName); // plot file name:
  129.  
  130.   // write input data for a check:
  131.   s_o << "-------------- these input data were read ---------------\n";
  132.   u->grid().print(s_o);  tip.print(s_o);
  133.   s_o << "Making plots at time points: "; time_points_for_plot.print (s_o);
  134.   s_o << "Flux function : " << flux_tp << "\ninitial function = ";
  135.   s_o << u0_tp << " uL=" << uL << '\n';
  136.   s_o << "---------------------------------------------------------\n";
  137. }
  138.  
  139. void ConsLaw:: setIC ()
  140. {
  141.   u_prev() = u0();
  142.   u_prev->valueIndex(grid->getBase(1)) = uL;
  143.  
  144.   // --- find the suitable dt ---
  145.   real dt;
  146.   real pmin, pmax;
  147.   p->minmax (pmin, pmax);
  148.  
  149.   if (flux_tp == "FluxLinear")
  150.     {
  151.       if (tip.Delta() > grid->Delta(1)/pmin)
  152.     dt = grid->Delta(1)/pmin;
  153.       else
  154.     dt = tip.Delta();
  155.     }
  156.   else if (flux_tp == "FluxSimpleBuckleyLeverett" || flux_tp == "FluxBurger")
  157.     dt = grid->Delta(1) / maxFluxDeriv (flux(),0,uL,0,0,0.1,0.1);
  158.   else if (flux_tp == "FluxBuckleyLeverett")
  159.     {
  160.       real pdelta;
  161.       if (pmin == pmax)
  162.     pdelta = 10;   // void infinite loop in maxFluxDeriv!
  163.       else if (pmax > pmin)
  164.     pdelta = pmax - pmin;
  165.       else
  166.     errorFP("ConsLaw::setIC","pmax=%f < pmin=%f ....",pmax,pmin);
  167.  
  168.       dt = grid->Delta(1) / maxFluxDeriv (flux(),0,uL,pmin,pmax,
  169.                       0.1,pdelta/10.0);
  170.     }
  171.   tip.setTimeStep (dt);    // store the computed dt in the TimePrm object
  172. }
  173.  
  174. void ConsLaw:: timeLoop ()
  175. {
  176.   tip.initTimeLoop();
  177.   while (!tip.finished())  {
  178.     tip.increaseTime();
  179.     solveAtThisTimeLevel ();
  180.     saveResults ();
  181.     updateDataStructures ();
  182.   }
  183. }
  184.  
  185. void ConsLaw:: solveAtThisTimeLevel ()
  186. {
  187.   scheme->update ();                    // execute the difference scheme
  188.   u->valueIndex(grid->getBase(1)) = uL; // set boundary conditions afterwards
  189. }
  190.  
  191. void ConsLaw:: saveResults ()
  192. {
  193.   if (time_points_for_plot.within (tip.getTime(), 0.4999*tip.Delta()))
  194.     DrawFD::makeCurvPlot (u(),file,"u(x) for t fixed","u",
  195.               oform("t=%g",tip.getTime()));
  196. }
  197.  
  198. void ConsLaw:: updateDataStructures () { u_prev() = u(); }
  199.